Objetivo: Fusionar productos de precipitación con datos de estaciones
La fusión de productos de precipitación con datos in-situ tiene como objetivo combinar y aprovechar la información proporcionada por ambos tipos de datos para mejorar la calidad y la resolución espacial de las estimaciones de precipitación.
Existen varias razones por las cuales se fusionan los productos de precipitación con los datos in-situ:
En el ejemplo de fusión de productos de precipitación con datos in-situ, se aplicará una metodología específica para combinar los datos y se mostrará cómo se pueden obtener estimaciones de precipitación mejoradas utilizando ambos tipos de datos.
La fusión de productos de precipitación con datos in-situ requiere considerar varios puntos:
homogeneidad de las fuentes de datos,
calidad y confiabilidad de los datos in-situ
Se pueden utilizar diferentes enfoques de fusión, como métodos estadísticos, interpolación espacial o modelos de fusión.
Baez-Villanueva et al., (2020)
En este ejemplo, fusionaremos dos productos de precipitación (CHIRPSv2 y PERSIANN-CDR) con estaciones in situ en Valparaíso, Chile.
Paso 1: Cargar los paquetes necesarios para trabajar con las estaciones:
El paquete RFmerge se puede instalar desde Github:
Paso 2: Cargar los datos de las estaciones y el shapefile del área de estudio. Para esto, almacenaremos las rutas de los archivos:
observations <- 'C:/Users/Data/Timeseries.zoo'
metadata <- 'C:/Users/Data/Metadata.csv'
shape <- 'C:/Users/Data/Valparaiso.shp'En este paso, estamos especificando las rutas de los archivos de datos de las estaciones y del shapefile del área de estudio. Asegúrate de proporcionar las rutas correctas a tus archivos correspondientes.
Podemos utilizar la función read.zoo del paquete zoo para importar las series temporales, la función read.csv para importar los metadatos y la función read_sf del paquete sf para importar el shapefile.
Paso 3: Ahora podemos cargar los datos diarios de los productos de precipitación y el modelo digital de elevación. Para esto, almacenaremos las rutas de los archivos correspondientes.
Podemos utilizar la función brick del paquete raster para importar los productos de precipitación (archivos multibanda) y la función raster para importar el DEM.
chirps <- brick(chirps)
persiann_cdr <- brick(persiann_cdr)
dem <- brick(dem)
print(chirps)
## class : RasterBrick
## dimensions : 40, 38, 1520, 243 (nrow, ncol, ncell, nlayers)
## resolution : 0.05, 0.05 (x, y)
## extent : -71.85, -69.95, -34, -32 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : CHIRPS5km.tif
## names : CHIRPS5km_1, CHIRPS5km_2, CHIRPS5km_3, CHIRPS5km_4, CHIRPS5km_5, CHIRPS5km_6, CHIRPS5km_7, CHIRPS5km_8, CHIRPS5km_9, CHIRPS5km_10, CHIRPS5km_11, CHIRPS5km_12, CHIRPS5km_13, CHIRPS5km_14, CHIRPS5km_15, ...
print(persiann_cdr)
## class : RasterBrick
## dimensions : 40, 38, 1520, 243 (nrow, ncol, ncell, nlayers)
## resolution : 0.05, 0.05 (x, y)
## extent : -71.85, -69.95, -34, -32 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : PERSIANNcdr5km.tif
## names : PERSIANNcdr5km_1, PERSIANNcdr5km_2, PERSIANNcdr5km_3, PERSIANNcdr5km_4, PERSIANNcdr5km_5, PERSIANNcdr5km_6, PERSIANNcdr5km_7, PERSIANNcdr5km_8, PERSIANNcdr5km_9, PERSIANNcdr5km_10, PERSIANNcdr5km_11, PERSIANNcdr5km_12, PERSIANNcdr5km_13, PERSIANNcdr5km_14, PERSIANNcdr5km_15, ...
print(dem)
## class : RasterBrick
## dimensions : 40, 38, 1520, 1 (nrow, ncol, ncell, nlayers)
## resolution : 0.05, 0.05 (x, y)
## extent : -71.85, -69.95, -34, -32 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : ValparaisoDEM5km.tif
## names : ValparaisoDEM5km
## min values : 11.4891
## max values : 5123.841Los dos archivos geotiff multibanda proporcionados no almacenan la fecha de la estimación de precipitación. Por lo tanto, antes de realizar cualquier análisis exploratorio, es necesario asignar nombres significativos a cada banda. Esto nos permitirá tener una mejor idea de las fechas asociadas a cada capa y facilitará el análisis posterior de los datos.
ldates <- hydroTSM::dip("1983-01-01", "1983-08-31")
## Please note that 'maptools' will be retired during October 2023,
## plan transition at your earliest convenience (see
## https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
## for guidance);some functionality will be moved to 'sp'.
## Checking rgeos availability: FALSE
names(chirps) <- ldates
names(persiann_cdr) <- ldatesA continuación, queremos visualizar las primeras y últimas seis filas de los metadatos espaciales:
head(metadata)
## Code lon lat NOM_REG NOM_PROV COD_COM NOM_COM
## 1 P5101005 -70.8000 -32.0836 Regin de Valparaso Los Andes 5304 San Esteban
## 2 P5111002 -71.0300 -32.1567 Regin de Valparaso Los Andes 5304 San Esteban
## 3 P5101006 -70.7833 -32.1836 Regin de Valparaso Los Andes 5304 San Esteban
## 4 P5100006 -70.7839 -32.2261 Regin de Valparaso Los Andes 5304 San Esteban
## 5 P5100005 -70.7100 -32.2286 Regin de Valparaso Los Andes 5304 San Esteban
## 6 P5111001 -71.1386 -32.2528 Regin de Valparaso Los Andes 5304 San Esteban
tail(metadata)
## Code lon lat NOM_REG NOM_PROV COD_COM NOM_COM
## 29 P5427006 -71.2144 -33.0986 Regin de Valparaso Los Andes 5304 San Esteban
## 30 P5510002 -71.5553 -33.1450 Regin de Valparaso Los Andes 5304 San Esteban
## 31 P5741002 -71.1467 -33.1686 Regin de Valparaso Los Andes 5304 San Esteban
## 32 P5530002 -71.6250 -33.5747 Regin de Valparaso Los Andes 5304 San Esteban
## 33 P5748003 -71.5106 -33.6344 Regin de Valparaso Los Andes 5304 San Esteban
## 34 P330030 -71.6142 -33.6550 Regin de Valparaso Los Andes 5304 San EstebanAhora podemos graficar la serie temporal diaria de precipitación para la primera estación (código: P5101005).
main <- paste("Precipitación diaria para la estación", metadata$Code[1])
ylab <- "Precipitación [mm]"
x_ts <- observations[,1]
#plot(x.ts, main=main, ylab= ylab, col="blue")
hydroTSM::hydroplot(x_ts, pfreq="o", main=main, ylab= ylab)
## [Note: pfreq='o' => ptype has been changed to 'ts']
grid()
Podemos graficar las estimaciones de precipitación acumulada para los primeros ocho meses de 1983 de CHIRPS y PERSIANN-CDR, y superponer los límites del área de estudio:
Paso 4: Convertir el objeto metadata en un objeto espacial. Para esto, lo transformatemos en un SpatialPointsDataFrame utilizando los campos de latitud y longitud, almacenados en las columnas lat y lon:
stations <- st_as_sf(metadata, coords = c('lon', 'lat'), crs = 4326)
print(stations)
## Simple feature collection with 34 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -71.625 ymin: -33.655 xmax: -70.4719 ymax: -32.0836
## Geodetic CRS: WGS 84
## First 10 features:
## Code NOM_REG NOM_PROV COD_COM NOM_COM
## 1 P5101005 Regin de Valparaso Los Andes 5304 San Esteban
## 2 P5111002 Regin de Valparaso Los Andes 5304 San Esteban
## 3 P5101006 Regin de Valparaso Los Andes 5304 San Esteban
## 4 P5100006 Regin de Valparaso Los Andes 5304 San Esteban
## 5 P5100005 Regin de Valparaso Los Andes 5304 San Esteban
## 6 P5111001 Regin de Valparaso Los Andes 5304 San Esteban
## 7 P5110003 Regin de Valparaso Los Andes 5304 San Esteban
## 8 P5111004 Regin de Valparaso Los Andes 5304 San Esteban
## 9 P5120004 Regin de Valparaso Los Andes 5304 San Esteban
## 10 P5200006 Regin de Valparaso Los Andes 5304 San Esteban
## geometry
## 1 POINT (-70.8 -32.0836)
## 2 POINT (-71.03 -32.1567)
## 3 POINT (-70.7833 -32.1836)
## 4 POINT (-70.7839 -32.2261)
## 5 POINT (-70.71 -32.2286)
## 6 POINT (-71.1386 -32.2528)
## 7 POINT (-70.9997 -32.2756)
## 8 POINT (-71.0792 -32.3058)
## 9 POINT (-71.2425 -32.3075)
## 10 POINT (-70.7528 -32.3408)Ahora podemos graficar el DEM con las ubicaciones de las estaciones pluviométricas que se utilizarán en el cálculo de RF-MEP, y superponiendo los límites del área de estudio y las estaciones:
Paso 5: Queremos producir un producto fusionado de precipitación desde el 1 de enero hasta el 31 de agosto de 1983 (243 días). Por lo tanto revisaremos la longitud de los productos.
Además, revisaremos la extensión espacial y la resolución espacial de los productos.
extent(chirps) == extent(persiann_cdr)
## [1] TRUE
res(chirps) == res(persiann_cdr)
## [1] TRUE TRUEPaso 6: Ahora vamos a fusionar los productos de precipitación, el DEM y los datos de las estaciones. Para esto, prepararemos una lista de covariables.
Si deseas que los archivos resultantes fusionados se escriban en disco, debes definir el directorio de salida (drty.out) antes de ejecutar RFmerge. Luego, puedes ejecutar la función RFmerge de la siguiente manera:
Si write2disk=TRUE, los siguientes archivos de salida se almacenarán en el directorio drty.out definido por el usuario:
Los objetos resultantes mencionados anteriormente se almacenarán dentro de drty.out de la siguiente manera:
Los objetos resultantes mencionados anteriormente se almacenarán dentro de drty.out de la siguiente manera:
Los objetos resultantes mencionados anteriormente se almacenarán dentro de drty.out de la siguiente manera:
rfmep_total <- sum(rfmep, na.rm = TRUE)
totals <- stack(chirps_total, persiann_total, rfmep_total)
names(totals) <- c('CHIRPSv2', 'PERSIANN_CDR', 'RF-MEP')
plot(totals)
Paso 7: utilizaremos el conjunto de datos de evaluación de observaciones de estaciones de medición de precipitación para evaluar el rendimiento de RF-MEP con observaciones que no se utilizaron en el procedimiento de fusión.
eval <- paste0(drty_out, "/Ground_based_data/Evaluation/Evaluation_ts.txt")
metadata_eval <- paste0(drty_out, "/Ground_based_data/Evaluation/Evaluation_metadata.txt")
eval <- read.zoo(eval, header = TRUE)
metadata_eval <- read.csv(metadata_eval)
stations_eval <- st_as_sf(metadata_eval, coords = c('lon', 'lat'), crs = 4326)Visualicemos una capa de RF-MEP junto con las estaciones de medición de precipitación y el área de estudio.
Primero, compararemos RF-MEP (y los dos productos de precipitación utilizados como covariables) con los datos de las estaciones de lluvia del conjunto de evaluación. Para extraer los valores de precipitación de RF-MEP en las ubicaciones de las estaciones, utilizaremos la función extract del paquete raster:
Ahora realizaremos una evaluación similar a la que hicimos en el Módulo 6 para evaluar el rendimiento de los productos de precipitación:
products <- list(chirps_ts, persiann_ts, rfmep_ts)
nproducts <- length(products)
nstations <- ncol(eval)
tmp <- rep(NA, nstations)
kges <- data.frame(ID=stations_eval[["Code"]], CHIRPS=tmp, PERSIANN_CDR=tmp, RF_MEP=tmp)
for (i in 1:nproducts) {
ldates <- time(eval)
lsim <- zoo(products[[i]], ldates)
kges[, (i+1)] <- hydroGOF::KGE(sim= lsim, obs= eval, method="2012")
}Ahora podemos visualizar los resultados:
print(kges)
## ID CHIRPS PERSIANN_CDR RF_MEP
## 1 P5120004 0.2456329 0.2568213 0.7087269
## 2 P5120003 0.4103445 0.3358265 0.8891823
## 3 P5414004 0.4492111 0.2953030 0.9106920
## 4 P5410007 0.4403373 0.1509242 0.9001426
## 5 P5423013 0.2988270 0.3147489 0.8865164
## 6 P5426004 0.3098025 0.3046178 0.5345708
## 7 P5427007 0.2280326 0.2647549 0.8525314Finalmente, se generará un boxplot que compare el desempeño, en términos de KGE’, del producto fusionado en comparación con los productos originales utilizados como covariables.
sres_cols <- c("powderblue", "palegoldenrod", "mediumseagreen")
boxplot(kges[,2:4], main = "Evaluación de KGE para enero - agosto de 1983",
xlab = "Productos de precipitación", ylab = "KGE'", ylim = c(0, 1), col=sres_cols)
grid()
legend("topleft", legend=c("CHIRPS", "PERSIANN-CDR", "RF-MEP"), col=sres_cols,
pch=15, cex=1.5, bty="n")
grid()